In [1]:
# necessary packages #

#using Pkg
#Pkg.add("Distances")
using Distributions
using Random
using Distances
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using ProgressMeter
using JLD2
In [2]:
include("../../../util2.j")
Out[2]:
colnorm (generic function with 1 method)
In [3]:
# unnecessary packages #

#using Pkg
#Pkg.add("UnicodePlots")
using UnicodePlots   # check the structure of the sparse matrix
using BenchmarkTools

using StatsPlots
using MCMCChains
using PrettyTables
In [4]:
#using Pkg
#Pkg.add("ProgressMeter");
In [5]:
@load "../data/sim3data.jld";
In [6]:
# priors #
p = size(X_ord)[2]; q = size(Y_ord)[2]; K = 2;
μβ = spzeros(p, q); inv_Vr = spzeros(p, p);
μΛ = spzeros(K, q); inv_VΛ = spzeros(K, K)
aΣ = 2; bΣ = fill(1.0, q);
inv_Lr = spzeros(p, p); inv_LΛ = spzeros(K, K);
aϕ = fill(2.0, K); bϕ = fill(33 / (2 * sqrt(2)), K);
In [11]:
# Some data preparations #

m = 10; n = length(S); 
# number of nearest neighbor                       
NN = BuildNN(coords_ord[S, :], m);                            # build nearest neighbor 
nnIndx_col = vcat(NN.nnIndx, 1:n);                            # the index of columns
nnIndx_row = zeros(Int64, 0);                                               
for i in 2:m
    nnIndx_row = vcat(nnIndx_row, fill(i, i-1))
end
nnIndx_row = vcat(nnIndx_row, repeat((m + 1):n, inner = m), 1:n);  # the index of rows

dim_invD = sum(index_S_M);
invD_yind = 1:dim_invD; invD_xind = 1:dim_invD;
Xtilde_indy_up = vcat([S .+ (ind - 1) * N for ind in 1:K]...);
nsam = length(perm_ind) + (K * n);
In [12]:
# preallocation #
μ_m = [Array{Float64, 2}(undef, length(M_ind[i]), q) for i in 1:q];
nIndx = length(NN.nnIndx);
A = [Array{Float64}(undef, nIndx) for i in 1:K];
D = [Array{Float64}(undef, n) for i in 1:K];
I_A = [spzeros(n, n) for i in 1:K];
A_new = [Array{Float64}(undef, nIndx) for i in 1:K];
D_new = [Array{Float64}(undef, n) for i in 1:K];
I_A_new = [spzeros(n, n) for i in 1:K];
Ystar = vcat(Y_ord[S, :], inv_Lr * μβ, inv_LΛ * μΛ);             # will be updated after imputing missing response
Xstar = vcat([X_ord[S, :] spzeros(n, K)], [inv_Lr spzeros(p, K)], 
    [spzeros(K, p) inv_LΛ]);      
bstar = fill(0.0, q); astar = aΣ + 0.5 * n;
μγstar = vcat(μβ, μΛ); invVγstar = fill(0.0, p + K, p + K);
Y_Xm = spzeros(n + p + K, q);
nsam = length(perm_ind) + (K * n);
Ytilde =  Array{Float64}(undef, nsam);
Xtilde = SparseMatrixCSC{Float64,Int64};
#precond_D = Array{Float64, 1}(undef, K * n);
lll = fill(1.0, (n, 1));

MCMC sampling algorithm Q1: priors for $\nu_i$ Q2: $\phi_i$ may not be consistant, since the order can change

In [13]:
# Preallocation for MCMC samples and Initalization #
N_sam = 20000;
N_pre_burn = Integer(trunc(0.5 * N_sam));
N_pre_adapt = Integer(trunc(0.25 * N_sam));
N_after_burn = N_sam - N_pre_burn;
ω_incp_sam = Array{Float32, 2}(undef, n, q);

ω_incp_sam_mean = fill(0.0, n, q);
ω_incp_sam_var = fill(0.0, n, q);
Y_m_sam_mean = [fill(0.0, length(M_ind[i])) for i in 1:length(M_ind)];
Y_m_sam_var = [fill(0.0, length(M_ind[i])) for i in 1:length(M_ind)];

F_sam = Array{Float64, 2}(undef, n, K);
Y_m_sam =  [Array{Float64, 1}(undef, length(M_ind[i])) for i in 1:q];
A_sam = Array{Float64, 2}(undef, N_sam, K); # acceptance rate
lh_old = 1; lh_new = 1;     # record the likelihood for updating ranges

ϕ_sam = Array{Float64, 2}(undef, K, N_sam + 1);

γ_sam = vcat(fill(0.0, p, q), fill(0.001, K, q));
γ_sam[(p + 1):(p + K), 1:K] = γ_sam[(p + 1):(p + K), 1:K] + I;
Σ_sam = fill(0.2, q);
ω_cov_sam = fill(0.0, q, q);
ϕ_sam[:, 1] = fill(6.0, K);

#precond_D = Array{Float64, 1}(undef, K * n);
inv_sqrt_Σ_diag = Array{Float64, 1}(undef, q);
RWM_scale = fill(0.1, K);
In [14]:
using DelimitedFiles
writedlm("../results/K2/γ_sam.csv", vcat(fill(0.0, 1, q), γ_sam), ", ");
writedlm("../results/K2/Σ_sam.csv", vcat(0.0, Σ_sam), ", ");
writedlm("../results/K2/ω_cov_sam.csv", [fill(0.0, q)], ", ");
writedlm("../results/K2/ϕ_sam.csv", vcat(0.0, ϕ_sam[:, 1]), ", ");
writedlm("../results/K2/A_sam.csv", 0.0, ", ");
In [15]:
# for loop for MCMC chain #
Random.seed!(123);
prog = Progress(N_sam, 1, "Computing initial pass...", 50)
for l in 1:N_sam
    # Build the matrix D_Sigma_o^{1/2} #
    inv_sqrt_Σ_diag = 1 ./ (sqrt.(Σ_sam));
    invD_ele = [x for ind in 1:q for x in fill(inv_sqrt_Σ_diag[ind], length(S_ind[ind]))];
    invD = sparse(invD_xind, invD_yind, invD_ele);
                
    if l == 1
        for k in 1:K
            getAD(coords_ord[S, :], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ϕ_sam[k, l], 0.5, 
                            A[k], D[k]);
            I_A[k] = sparse(nnIndx_row, nnIndx_col, vcat(-A[k], ones(n)));
        end
    end
    Ytilde = vcat(invD * vcat([Y_ord[S_ind[ind], ind] - X_ord[S_ind[ind], :] * 
                            γ_sam[1:p, ind] for ind in 1:q]...), zeros(K * n));
    Xtilde = vcat(invD * kron(sparse(transpose(γ_sam[(p + 1):(p + K), :])), 
                            sparse(1:N, 1:N, ones(N)))[obs_ind, Xtilde_indy_up],
             blockdiag([Diagonal(1 ./ sqrt.(D[ind])) * I_A[ind] for ind in 1:K]...));
                  
    # use LSMR to generate sample of F #       
    #Precond_D = colnorm(Xtilde);             
    #F_sam = reshape(Diagonal(1 ./ Precond_D) * lsmr(Xtilde * Diagonal(1 ./ Precond_D), 
    #                     collect(Ytilde) + rand(Normal(), nsam)), :, K);                
    F_sam = reshape(lsmr(Xtilde, collect(Ytilde) + rand(Normal(), nsam)), :, K);
                
    # impute missing response  over S#
    Xstar[1:n, (p + 1):(p + K)] = F_sam;        # update matrix Xstar with F
    if(l > N_pre_burn) # only save ω_incp_sam after burn-in
        ω_incp_sam = F_sam * γ_sam[(p + 1):(p + K), :] + lll * transpose(γ_sam[1, :]); 
        ω_incp_sam_mean = ω_incp_sam_mean + (ω_incp_sam ./ N_after_burn);
        ω_incp_sam_var = ω_incp_sam_var + ((ω_incp_sam.^2) ./ N_after_burn);  
        ω_cov_sam = cov(ω_incp_sam);
    else
        ω_cov_sam = cov(F_sam * γ_sam[(p + 1):(p + K), :]);
    end                    
    io1 = open("../results/K2/ω_cov_sam.csv", "a" ); # covariance of latent process
    writedlm(io1, ω_cov_sam, ", ");
    close(io1);            
           
    # impute missing response  over S#            
    for ind in 1:q
        Y_m_sam[ind] = Xstar[M_Sind[ind], :] * γ_sam[:, ind]+   
            rand(Normal(0.0, sqrt(Σ_sam[ind])), length(M_ind[ind]));
    end
         
    if (l > N_pre_burn)  # only save imputed Y after burn-in
        for ind in 1:q
            Y_m_sam_mean[ind] = Y_m_sam_mean[ind] + (Y_m_sam[ind] ./ N_after_burn);
            Y_m_sam_var[ind] = Y_m_sam_var[ind] + ((Y_m_sam[ind].^2) ./ N_after_burn);
        end
    end
                
                        
    # use MNIW to sample γ Σ #
    for ind in 1:q
        Ystar[M_Sind[ind], ind] = Y_m_sam[ind]   # update Ystar with imputed response
    end
    
    invVγstar = cholesky(Xstar'Xstar);
    mul!(μγstar, transpose(Xstar), Ystar); μγstar = invVγstar.U \ (invVγstar.L \ μγstar);
    Y_Xm = Ystar - Xstar * μγstar;      # maybe improve?
    bstar = [bΣ[ind] + 0.5 * (norm(Y_Xm[:, ind])^2) for ind in 1:q];
    Σ_sam = [rand(InverseGamma(astar, bstar[ind]), 1)[1] for ind in 1:q];          # sample Σ
    γ_sam = (invVγstar.U \ reshape(rand(Normal(), (p + K) * q), (p + K), q)) * 
                    Diagonal(sqrt.(Σ_sam)) + μγstar;          # sample γ    
    io2 = open("../results/K2/Σ_sam.csv", "a" );
    writedlm(io2, Σ_sam, ", ");
    close(io2); 
    io3 = open("../results/K2/γ_sam.csv", "a" );
    writedlm(io3, γ_sam, ", ");
    close(io3)
                
    # use adaptive metropolis-hasting to update range
    if l > 3 && l < N_pre_adapt
        RWM_scale = [sqrt(2.38^2 * var(ϕ_sam[i, Integer(floor(l / 2)):l], 
                                corrected=true) * 0.95^2 + 0.05^2 * 0.1^2) for i in 1:K];
    end
    # use metropolis-hasting to update range
    ϕ_sam[:, l + 1] = [ϕ_sam[i, l] + RWM_scale[i] * rand(Normal(), 1)[1] for i in 1:K]; # propose next sample point
  
    for i in 1:K
        if ϕ_sam[i, l + 1] > 0.0
            lh_old = -0.5 * (sum(log.(D[i])) + norm((I_A[i] * F_sam[:, i]) ./ sqrt.(D[i]))^2) + 
               loglikelihood(Gamma(aϕ[i], bϕ[i]), [ϕ_sam[i, l]]);
            getAD(coords_ord[S, :], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ϕ_sam[i, l + 1], 
                0.5, A_new[i], D_new[i]);
            I_A_new[i] = sparse(nnIndx_row, nnIndx_col, vcat(-A_new[i], ones(n)));
            lh_new = -0.5 * (sum(log.(D_new[i]))  + norm((I_A_new[i] * F_sam[:, i]) ./ sqrt.(D_new[i]))^2) +
                   loglikelihood(Gamma(aϕ[i], bϕ[i]), [ϕ_sam[i, l + 1]]);     
            A_sam[l, i] = min(exp(lh_new - lh_old), 1.0);
            if rand(1)[1] < A_sam[l, i]
                I_A[i] = copy(I_A_new[i]); D[i] = copy(D_new[i]);        # update and update the corresponding I_A D
            else 
                ϕ_sam[i, l + 1] = ϕ_sam[i, l];
            end
        else 
            A_sam[l, i] = 0.0;
            ϕ_sam[:, l + 1] = ϕ_sam[:, l];   
        end
    end                   
    
    io4 = open("../results/K2/ϕ_sam.csv", "a" );
    writedlm(io4, ϕ_sam[:, l + 1], ", ");
    close(io4); 
    io5 = open("../results/K2/A_sam.csv", "a" );
    writedlm(io5, A_sam[l, :], ", ");
    close(io5)
                
    next!(prog) # monitor the progress
end
ω_incp_sam_var = (ω_incp_sam_var - ω_incp_sam_mean.^2) * (N_after_burn / (N_after_burn - 1));
Y_m_sam_var = [(Y_m_sam_var[ind] - Y_m_sam_mean[ind].^2) * 
               (N_after_burn / (N_after_burn - 1)) for ind in 1:q];

Computing initial pass...100%|██████████████████████████████████████████████████| Time: 0:07:02:07
In [16]:
round.([mean(A_sam[(N_pre_burn + 1):N_sam, i]) for i in 1:K], digits = 5)
Out[16]:
2-element Array{Float64,1}:
 0.02378
 0.33766
In [17]:
RWM_scale
Out[17]:
2-element Array{Float64,1}:
 1.8056616534482128 
 0.34536277973900764

Posterior prediction

In [18]:
# prediction preparison
#N_pre_burn = Integer(trunc(0.75 * N_sam));
#M_ind = setdiff(1:N, S); NM = length(M_ind)
#F_M_sam = Array{Float64, 3}(undef, NM, K, N_sam - N_pre_burn + 1);
#Y_M_sam = Array{Float64, 3}(undef, NM, q, N_sam - N_pre_burn + 1);

# construct Atilde Dtilde #

#using RCall
#@rput coords_ord
#@rput S
#@rput m
#R"""
#library("RANN")
#nn_mod_ho <- nn2(t(coords_ord[, S]), t(coords_ord[, -S]), k = m)
#"""
#@rget nn_mod_ho
#Atilde = Array{Float64}(undef, NM * m); Dtilde = Array{Float64}(undef, NM);
#MnnIndxLU = collect(1:m:(NM * m + 1));
#MnnIndx = vec(nn_mod_ho[:nn_idx]');
#Mnndists = vec(nn_mod_ho[:nn_dists]')
In [19]:
#for i in N_pre_burn:N_sam
#    for j in 1:K
#        # update F
#        getAD(coords_ord[:, S], MnnIndx, vec(nn_mod_ho[:nn_dists]'), MnnIndxLU, 
#             ϕ_sam[j, i + 1], 0.5, Atilde, Dtilde)
#        AtildeM = sparse(repeat(1:NM, inner = m), MnnIndx, Atilde, NM, n);
#        F_M_sam[:, j, (i - N_pre_burn + 1)] = AtildeM * F_sam[:, j, i] + sqrt.(Dtilde) .* rand(Normal(), NM)
#    end 
#    # update Y
#    Y_M_sam[:, :, (i - N_pre_burn + 1)] = X_ord[M_ind, :] * γ_sam[1:p, :, i + 1] + 
#        F_M_sam[:, :, (i - N_pre_burn + 1)] *  γ_sam[(p + 1):(p + K), :, i + 1] + 
#        transpose(rand(MvNormal(Σ_sam[:, :, i + 1]), NM))
#end
In [20]:
#load data
using CSV
γ_sam = convert(Matrix{Float64}, CSV.read("../results/K2/γ_sam.csv"));
ind_γ_sam = 1: (p + K) :((p + K) * N_sam + 1);
Σ_sam = convert(Matrix{Float64}, CSV.read("../results/K2/Σ_sam.csv"));
ind_Σ_sam = 1: q :(q * N_sam + 1);
ω_cov_sam = convert(Matrix{Float64}, CSV.read("../results/K2/ω_cov_sam.csv"));
ind_ω_cov_sam = 1: q :(q * (N_sam - 1) + 1);

MCMC Chain check

In [21]:
β_pos_sam = Array{Float64, 3}(undef, N_sam + 1, (p -1) * q, 1);
β_pos_sam[:, :, 1] = hcat([γ_sam[ind_γ_sam .+ (i - 1), j] for i in 2:p, j in 1:q]...);
β_chain = Chains(β_pos_sam);
pβ = plot(β_chain)
Out[21]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -6 -5 -4 -3 -2 -1 0 Param1 Iteration Sample value -6 -5 -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 8 Param2 Iteration Sample value 0 2 4 6 8 0.0 0.5 1.0 1.5 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 Param3 Iteration Sample value 0 1 2 3 4 0.0 0.2 0.4 0.6 0.8 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 8 Param4 Iteration Sample value 0 2 4 6 8 0.0 0.2 0.4 0.6 0.8 Param4 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 Param5 Iteration Sample value 0 1 2 3 4 0.0 0.3 0.6 0.9 1.2 Param5 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -12 -10 -8 -6 -4 -2 0 Param6 Iteration Sample value -12 -10 -8 -6 -4 -2 0 0.0 0.3 0.6 0.9 1.2 Param6 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -3 -2 -1 0 1 Param7 Iteration Sample value -4 -3 -2 -1 0 1 0.0 0.2 0.4 0.6 Param7 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -2 -1 0 1 2 Param8 Iteration Sample value -3 -2 -1 0 1 2 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Param8 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -8 -6 -4 -2 0 Param9 Iteration Sample value -8 -6 -4 -2 0 0.0 0.2 0.4 0.6 Param9 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -6 -5 -4 -3 -2 -1 0 Param10 Iteration Sample value -6 -5 -4 -3 -2 -1 0 0.0 0.2 0.4 0.6 Param10 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 5 6 Param11 Iteration Sample value 0 1 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 Param11 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0.0 2.5 5.0 7.5 10.0 Param12 Iteration Sample value 0.0 2.5 5.0 7.5 10.0 0.0 0.2 0.4 0.6 0.8 Param12 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 5 6 Param13 Iteration Sample value 0 1 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 Param13 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -10.0 -7.5 -5.0 -2.5 0.0 Param14 Iteration Sample value -10.0 -7.5 -5.0 -2.5 0.0 0.0 0.2 0.4 0.6 Param14 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -5 -4 -3 -2 -1 0 Param15 Iteration Sample value -5 -4 -3 -2 -1 0 0.0 0.2 0.4 0.6 Param15 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 Param16 Iteration Sample value 0 2 4 6 0.0 0.2 0.4 0.6 Param16 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 8 Param17 Iteration Sample value 0 2 4 6 8 0.0 0.2 0.4 0.6 Param17 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 Param18 Iteration Sample value 0 2 4 6 8 0.0 0.2 0.4 0.6 Param18 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -5 -4 -3 -2 -1 0 Param19 Iteration Sample value -6 -5 -4 -3 -2 -1 0 0.0 0.2 0.4 0.6 0.8 Param19 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -8 -6 -4 -2 0 Param20 Iteration Sample value -8 -6 -4 -2 0 0.0 0.2 0.4 0.6 0.8 Param20 Sample value Density
In [22]:
β
Out[22]:
3×10 Array{Float64,2}:
  1.0  -1.0    1.0  -0.5   2.0  -1.5   0.5   0.3  -2.0   1.5
 -5.0   2.0    3.0  -2.0  -6.0   4.0   5.0  -3.0   6.0  -4.0
  8.0   6.9  -12.0   0.0  -4.0   7.7  -8.8   3.3   6.6  -5.5
In [23]:
Λ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, K * q, 1);
Λ_pos_sam[:, :, 1] = hcat([γ_sam[ind_γ_sam .+ (p + i - 1), j] for i in 1:K, j in 1:q]...)
Λ_chain = Chains(Λ_pos_sam);
#pΛ = plot(Λ_chain)
In [24]:
Λ;
In [25]:
ϕ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, K, 1);
ϕ_pos_sam[:, :, 1] = hcat([ϕ_sam[i, :] for i in 1:K]...);
ϕ_chain = Chains(ϕ_pos_sam);
pϕ = plot(ϕ_chain)
Out[25]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 4 6 8 10 Param1 Iteration Sample value 0.0 2.5 5.0 7.5 10.0 0.0 0.2 0.4 0.6 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 5 6 Param2 Iteration Sample value 0 1 2 3 4 5 6 0 1 2 3 Param2 Sample value Density
In [26]:
Σ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, q, 1);
Σ_pos_sam[:, :, 1] = hcat([Σ_sam[ind_Σ_sam .+ (ind - 1)] for ind in 1:q]...);
Σ_chain = Chains(Σ_pos_sam);
pΣ = plot(Σ_chain)
Out[26]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 8 Param1 Iteration Sample value 0 2 4 6 8 0.0 0.2 0.4 0.6 0.8 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 Param2 Iteration Sample value 0 5 10 15 0.0 0.1 0.2 0.3 0.4 0.5 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 10 20 30 Param3 Iteration Sample value 0 10 20 30 0.0 0.2 0.4 0.6 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 10 20 30 Param4 Iteration Sample value 0 10 20 30 0.0 0.1 0.2 0.3 Param4 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 25 Param5 Iteration Sample value 0 5 10 15 20 25 0.0 0.1 0.2 0.3 Param5 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 25 Param6 Iteration Sample value 0 5 10 15 20 25 0.0 0.1 0.2 0.3 0.4 Param6 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 25 Param7 Iteration Sample value 0 5 10 15 20 25 0.0 0.1 0.2 0.3 Param7 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 10 20 30 Param8 Iteration Sample value 0 10 20 30 0.0 0.1 0.2 0.3 Param8 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 25 Param9 Iteration Sample value 0 10 20 30 0.0 0.1 0.2 0.3 Param9 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 Param10 Iteration Sample value 0 5 10 15 20 0.0 0.1 0.2 0.3 0.4 Param10 Sample value Density
In [27]:
sqrt_Σ_diag.^2
Out[27]:
10-element Array{Float64,1}:
 0.5000000000000001 
 1.0                
 0.4                
 2.0000000000000004 
 0.29999999999999993
 2.5000000000000004 
 3.5                
 0.45               
 1.4999999999999998 
 0.5000000000000001
In [28]:
ω_cov_pos_sam = Array{Float64, 3}(undef, N_sam, 2 * 2, 1);
ω_cov_pos_sam[:, :, 1] = hcat(ω_cov_sam[ind_ω_cov_sam, 1], ω_cov_sam[ind_ω_cov_sam, 2], 
    ω_cov_sam[ind_ω_cov_sam .+ 1, 1], ω_cov_sam[ind_ω_cov_sam .+ 1, 2]);
ω_cov_chain = Chains(ω_cov_pos_sam);
pωcov = plot(ω_cov_chain)
Out[28]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 14 16 18 20 22 24 26 Param1 Iteration Sample value 15 18 21 24 27 0.0 0.1 0.2 0.3 0.4 0.5 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 3 4 5 6 7 Param2 Iteration Sample value 2 3 4 5 6 7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 3 4 5 6 7 Param3 Iteration Sample value 2 3 4 5 6 7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2.5 5.0 7.5 10.0 12.5 Param4 Iteration Sample value 2.5 5.0 7.5 10.0 12.5 0.00 0.25 0.50 0.75 1.00 Param4 Sample value Density
In [29]:
covω = cov(ω_ord[S, :])
Out[29]:
10×10 Array{Float64,2}:
  27.3056    5.61976     6.67839   …   3.12784   -13.7544    -13.057  
   5.61976  17.1003      7.39449      -0.188053   -6.98366    -7.9723 
   6.67839   7.39449    28.0238        1.15105   -15.5356    -14.9628 
   1.84707   1.93173   -10.7337       -3.4337      0.41701     4.35427
  -5.90397   4.82641     0.199311      1.52289    -0.756625    1.53052
  16.4464    5.06729    17.1642    …   6.00069   -12.0785    -23.633  
 -14.2019   -3.08495     1.25         -3.33428     1.39374     6.75592
   3.12784  -0.188053    1.15105      26.2392      9.25874    -4.51342
 -13.7544   -6.98366   -15.5356        9.25874    33.0245     10.3836 
 -13.057    -7.9723    -14.9628       -4.51342    10.3836     30.874
In [30]:
corω = cor(ω_ord[S, :])
Out[30]:
10×10 Array{Float64,2}:
  1.0         0.26007      0.241425    …   0.116854    -0.458035   -0.449699
  0.26007     1.0          0.337787       -0.00887776  -0.293875   -0.346964
  0.241425    0.337787     1.0             0.0424478   -0.510677   -0.508691
  0.0622972   0.082329    -0.357353       -0.11814      0.012789    0.138111
 -0.2279      0.235423     0.00759441      0.059968    -0.0265576   0.055561
  0.52487     0.204353     0.540712    …   0.195359    -0.350511   -0.709299
 -0.525134   -0.144143     0.0456241      -0.125769     0.046861    0.234929
  0.116854   -0.00887776   0.0424478       1.0          0.314528   -0.158575
 -0.458035   -0.293875    -0.510677        0.314528     1.0         0.325187
 -0.449699   -0.346964    -0.508691       -0.158575     0.325187    1.0

Posterior Inference

In [31]:
ω_cov_pos_sam_mean = [mean(ω_cov_sam[ind_ω_cov_sam .+ (i - 1), j][(N_pre_burn + 1):N_sam]) 
    for i in 1:q, j in 1:q]
Out[31]:
10×10 Array{Float64,2}:
  24.7402    4.58562     5.75077   …   3.68037   -12.2914   -14.1325 
   4.58562   2.3776      7.11216       0.902887   -5.18498   -5.9343 
   5.75077   7.11216    25.8483        1.74753   -14.6316   -16.6989 
   1.78399  -2.50117   -11.0635       -0.154051    4.62789    5.26954
  -6.25421  -0.849976   -0.208192     -0.885428    2.50776    2.89182
  18.1603    6.91686    18.5954    …   3.22932   -15.9315   -18.2655 
 -15.299    -1.55304     1.65012      -2.08946     5.09941    5.89802
   3.68037   0.902887    1.74753       0.637972   -2.2537    -2.59508
 -12.2914   -5.18498   -14.6316       -2.2537     11.8193    13.4693 
 -14.1325   -5.9343    -16.6989       -2.59508    13.4693    15.4716
In [32]:
ω_corr_sam = [(Diagonal([1 / sqrt(ω_cov_sam[ind_ω_cov_sam[l] .+ (i - 1), i]) for i in 1:q]) * 
    ω_cov_sam[ind_ω_cov_sam[l] .+ (1:q) .- 1, 1:q] * 
        Diagonal([1 / sqrt(ω_cov_sam[ind_ω_cov_sam[l] .+ (i - 1), i]) for i in 1:q])) for l in 1:N_sam];
In [33]:
ω_corr_sam_mean = [mean([ω_corr_sam[i][j , k] for i in (N_pre_burn + 1):N_sam]) 
    for j in 1:q, k in 1:q]
Out[33]:
10×10 Array{Float64,2}:
  1.0        0.600319   0.227446   …   0.945777   -0.719612  -0.723158
  0.600319   1.0        0.910542       0.752114   -0.982662  -0.982427
  0.227446   0.910542   1.0            0.4409     -0.838121  -0.835854
  0.151742  -0.691998  -0.924885      -0.0849732   0.573171   0.569653
 -0.972842  -0.428959  -0.031541      -0.875038    0.56514    0.569217
  0.782191   0.96413    0.783587   …   0.884469   -0.993655  -0.995016
 -0.945125  -0.310721   0.0995629     -0.819224    0.456446   0.460916
  0.945777   0.752114   0.4409         1.0        -0.840488  -0.844377
 -0.719612  -0.982662  -0.838121      -0.840488    1.0        0.997897
 -0.723158  -0.982427  -0.835854      -0.844377    0.997897   1.0
In [34]:
# check the plot of the data 
using RCall
@rput ω_corr_sam_mean
R"""
library(corrplot)
library(corrgram)

colnames(ω_corr_sam_mean) <- c(1:10)
rownames(ω_corr_sam_mean) <- c(1:10)

width <- 720
height <- 720
pointsize <- 16


png(paste("../../pics/sim3_factor_corr_K2_plot.png", sep = ""), 
    width = width, height = height, pointsize = pointsize, family = "Courier")
par(mfrow = c(1, 1))
#corrplot(ω_corr_sam_mean, method="number", type = "upper", diag = FALSE, 
#         addshade = "negative", number.digits = 2, tl.pos = "td")

corrgram(ω_corr_sam_mean, order=FALSE, lower.panel=panel.shade, gap = 0.2,
         upper.panel=panel.pie, text.panel=panel.txt, main=" ",
         col.regions = colorRampPalette(c( "darkseagreen3",
                                           "white", "cadetblue3")))
dev.off()
"""

┌ Warning: Package RCall does not have AxisArrays in its dependencies:
│ - If you have RCall checked out for development and have
│   added AxisArrays as a dependency but haven't updated your primary
│   environment's manifest file, try `Pkg.resolve()`.
│ - Otherwise you may need to report an issue with RCall
│ Loading AxisArrays into RCall from project dependency, future warnings for RCall are suppressed.
└ @ nothing nothing:840
┌ Warning: RCall.jl: corrplot 0.84 loaded
└ @ RCall /home/lu/.julia/packages/RCall/ffM0W/src/io.jl:113
┌ Warning: RCall.jl: Registered S3 method overwritten by 'seriation':
│   method         from 
│   reorder.hclust gclus
└ @ RCall /home/lu/.julia/packages/RCall/ffM0W/src/io.jl:113
Out[34]:
RObject{IntSxp}
null device 
          1
In [35]:
# CVG-slope #
count_slope = 0.0
for i in 2:p
    for j in 1:q
        if ((quantile(γ_sam[ind_γ_sam .+ (i - 1), j][(N_pre_burn + 1):(N_sam + 1)], [0.025])[1] < 
                β[i, j] ) && (quantile(
                        γ_sam[ind_γ_sam .+ (i - 1), j][(N_pre_burn + 1):(N_sam + 1)], [0.975])[1] > 
                β[i, j] ))
            count_slope = count_slope + 1.0;
        end
        
    end
end
count_slope
Out[35]:
18.0
In [36]:
# CVL #
count_ω_incp = fill(0.0, q);
for j in 1:q
    for i in 1:n
        count_ω_incp[j] = count_ω_incp[j] + 
        (((ω_incp_sam_mean[i, j] - 1.96 * sqrt(ω_incp_sam_var[i, j])) < ω_incp_obs[S[i], j]) && 
            ((ω_incp_sam_mean[i, j] + 1.96 * sqrt(ω_incp_sam_var[i, j])) > ω_incp_obs[S[i], j]))
    end
end
count_ω_incp;
In [37]:
print(round.(count_ω_incp ./ n, digits = 5))

[0.9175, 0.215, 0.935, 0.2675, 0.18917, 0.43167, 0.36333, 0.1575, 0.33083, 0.3675]
In [38]:
sum(count_ω_incp) / (q*n)
Out[38]:
0.4175
In [39]:
# CVG #
count_Y_M = fill(0.0, q);
Y_m_pos_qt = [Array{Float64, 2}(undef, length(M_ind[ind]), 3) for ind in 1:q];

for i in 1:q
    for j in 1:length(M_ind[i])
        count_Y_M[i] = count_Y_M[i] + (((Y_m_sam_mean[i][j] - 
                1.96 * sqrt(Y_m_sam_var[i][j])) < Y_ord[M_ind[i][j], i]) && 
         ((Y_m_sam_mean[i][j] + 1.96 * sqrt(Y_m_sam_var[i][j])) > 
                Y_ord[M_ind[i][j], i]))
    end
end
print(round.(count_Y_M ./ 200, digits = 4));
print(round.(mean(count_Y_M ./ 200), digits = 4))

[0.925, 0.94, 0.93, 0.945, 0.945, 0.925, 0.945, 0.945, 0.965, 0.955]0.942
In [40]:
# RMSPE #
MSPE = (sum([sum((Y_m_sam_mean[i] - Y_ord[S[M_Sind[i]], i]).^2) for i in 1:q])) / 
    (sum([length(M_Sind[i]) for i in 1:q]))
RMSPE = sqrt(MSPE); 
print(round.([sqrt((sum((Y_m_sam_mean[i] - Y_ord[S[M_Sind[i]], i]).^2)) / 
        length(M_Sind[i])) for i in 1:q], digits = 4));
print(round(RMSPE, digits = 4));

[2.5998, 4.0007, 2.7377, 5.5204, 4.8977, 4.5303, 4.7247, 5.0411, 4.7013, 4.1487]4.3852
In [42]:
@save "../results/K2/Factor_mean_var_K2.jld" ω_incp_sam_mean ω_incp_sam_var Y_m_sam_mean Y_m_sam_var Y_ord S_ind M_Sind K p q N_sam N_pre_burn